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ABSTRACT 

We test statistically the hypothesis that radio pulsar glitches result from an 
avalanche process, in which angular momentum is transferred erratically from 
the flywheel-like superfluid in the star to the slowly decelerating, solid crust 
via spatially connected chains of local, impulsive, threshold-activated events, so 
that the system fluctuates around a self-organised critical state. Analysis of the 
glitch population (currently 285 events from 101 pulsars) demonstrates that the 
size distribution in individual pulsars is consistent with being scale invariant, as 
expected for an avalanche process. The measured power-law exponents fall in 
the range —0.13 < a < 2.4, with a ^ 1.2 for the youngest pulsars. The waiting- 
time distribution is consistent with being exponential in seven out of nine pulsars 
where it can be measured reliably, after adjusting for observational limits on the 
minimum waiting time, as for a constant-rate Poisson process. PSR J0537— 6910 
and PSR J0835— 4510 are the exceptions; their waiting-time distributions show 
evidence of quasiperiodicity. In each object, stationarity requires that the rate 
A equals — ez>/(Ai/), where z> is the angular acceleration of the crust, (Az/) is 
the mean glitch size, and ez> is the relative angular acceleration of the crust and 
superfluid. Measurements yield e < 7 x 10~^ for PSR J0358-I-5413 and e < 1 
(trivially) for the other eight objects, which have a < 2. There is no evidence 
that A changes monotonically with spin-down age. The rate distribution itself is 
fitted reasonably well by an exponential for A > 0.25 yr^^, with (A) = 1.3lQ'g yr~^. 
For A < 0.25 yr^^, its exact form is unknown; the exponential overestimates the 
number of glitching pulsars observed at low A, where the limited total observation 
time exercises a selection bias. In order to reproduce the aggregate waiting-time 
distribution of the glitch population as a whole, the fraction of pulsars with 
A > 0.25 yr^^ must exceed ~ 70 per cent. 
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Subject headings: dense matter — pulsars: general — stars: interior — stars: 
neutron — stars: rotation 



Introduction 



Glitches are tiny, impulsive, randomly timed increases in the spin frequency of a 
rotation-powered pulsar, sometimes accompanied by an impulsive change in the frequency 
derivative z>. They are to be distinguished from timing noise, a type of rotational irregularity 
where pulse arrival times wander continuously, although there is evi dence that timing noise 



is the cumulative result o f frequent microglitches in certain pulsars (ICordes &: Downs! Il985 
D'Alessandro et alj[l995r i. 



At the time of writing, 285 glitches in total have been detected in 101 objects (~ 6% 
of the known radio pulsar population), the majority in the last four years, facilitated by 
the Parkes Multibeam S urvey, refined multifrequency ephe r nerides, and better interfer- 
ence rejection algorithms (IHobbsl l2002l : iKrawczyk et al.l l2003t iKramer fc Lynd l2005l : iLewis 



20051 : IJanssen fc Stappers 



20061 ). Efforts to analyse the data statis t ically have focused on 



the correlation of glit c h activity with age (jMcKenna fc Lvnd Il990l : IShemar fc Lyne 



2006 



200ll ) , the size distribution (IMorley fc Garcia-Pelavo 



tion between glitch si zes and waiting times (IWang et al. 



1993a 



2000 



Per alt a 



Wong et al 



1996 



Urama fc Ok^ Il999l: iLvne et all boool : IWang et al.lboool ) and Re ynolds number JPeraltal 



Melatos fc PeraltalboOTh. the post-glitch relaxation tim e -scale JWang et al.|[200ol : IWong et~al 



2006), and the correla- 



2001: Middleditch et al. 



20061 : |Peraltall2006l ). iHobbd (120021 ) reviewed the role of observational selection effects. 



Most glitching pulsars (65%) have been seen to glitch once, but a minority glitch re- 
peatedly; the current record holder is PSR J1740— 3015, with 33 glitches. Of those objects 
which glitch repeatedly, most do so at unpredictable intervals, but two (PSR J0537— 6910 
and Vela) are qua siperiodic; Vela, in particular, has been likened to a relaxation oscillator 
( iLyne et al.lll996l ). The fractional increase in u spans seven decades (3 x 10~^^ < Au < 



2 X 10 ^) across the glitch population and as many as four decades in a single object (e.g. 
7 X 10-1° < Az/ < 2 X 10-6 in PSR J1740-3015). The spin-down age = -z//(2z>) of glitch- 
ing pulsars spans four decades, from 1 x 10^ yr to 3 x 10'' yr. In many respects, therefore, the 
glitch phenomenon is scale invariant. This striking property invites physical interpretation. 

Theories of pulsar glitches have focused mainly on the local microphysics of the super- 
fluid in the ste llar interior and its coupling to the solid crust, for example the strength of 
vortex pinning (lAnderson fc Itohlll975l : |Joneslll998l ) , the rate of vortex creep (ILink fc Epstein 
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1996 ) , or the conditions for ex citing; superfluid turbulence ( Peralta et all2005 , 2006 : Melatos fc Peralta 
20071 : lAndersson et al.l 120071 ). Ultimately, however, the local microphysics must be synthe- 
sized with the global, collective dynamics in order to make full contact with observational 
data. (Likewise, a practical model of earthquakes must synthesize the microphysics of rock 
fracture with the macrodynamics of interacting tectonic plates.) For example, if approxi- 
mately 10^^(Az//l Hz) vortices unpin from crustal lattice sites in sympathy during a glitch, 
they must communicate rapidly across distances much greater than their separation. How? 
And why does the number that unpin fluctuate so dramatically (by up to four orders of mag- 
nitude) from glitch to glitch in a single pulsar, while always amounting to a small fraction 
(Az//z/) of the total? 

Such collective, scale invariant behavior is a generic feature of a class of natural and syn- 
thetic far-from-equilibrium systems, called self-organized critical systems, that are discrete, 
interaction dominated, and slowly driven, and that adjust internally via erratic, spatially con - 
nected avalanches of local, impulsive, threshold-activated, relaxation events (|Jensenlll998l ). 
Such systems fluctuate around a stationary state towards which they evolve spontaneously, in 
which global driving balances local relaxation on average over the long term. The archetype 
of a self-organized critical system is the sandpile (IBak et al.l 119871 ). 



In this paper, we study pu lsar glitches as an avalanche process, as first proposed by 
Morley fc Garcfa-Pelayd (Il993al ). After reviewing self-organised criticality in ^ we define 
the statistical sample on which our study is based (§!]) and analyze the observed distribu- 
tion of glitch sizes (§!]) and waiting times (^. Some implications for glitch physics are 
explored in ^ We only include radio pulsars in the sample, to preserve its homogeneity, 
even though glitches have now been observed in anomalous X-ray pulsars (magnetars) as 
well dPairOsso et~aDl2003l : E<aspi fc Gavriilll2003h . 



2. Avalanche dynamics 



A sy stern in a self-organised critical state exhibits the following distinguishing features 
Jjensenlll998l ). 



1. It is composed of many discrete, mutually interacting elements, whose motions are 
dominated by local (e.g. nearest neighbor) rather than global (e.g. mean field) forces. 

a 



^Tectonic plates, or grains of sand in a pile, are terrestrial examples of interacting elements. 



-4- 



Each element moves when the local force exceeds a threshold (stick-slip motion). Hence 
stress accumulates sustainedly at certain random locations while relaxing quickly else- 
where; at any instant, the system houses numerous metastable stress reservoirs, sepa- 
rated by relaxed zones. 

An external force drives the system slowly, in the sense that elements adjust to local 
forces rapidly compared to the driver time-scale. Combined with local thresholds, this 
ensures that the system evolves quasistatically through a history-dependent sequence 
of metastable states (a huge number of which are available). 

Transitions from one metastable state to the next occur via avalanches: spatially con- 
nected chains of local equilibration events, in which one element relaxes and redis- 
tributes some local stress to its neighbors, which in turn can exceed their thresholds 
and relax (knock-on effect). The duration of even the largest avalanches is short com- 
pared to the driving time-scale (see previous point). 

Avalanches have no preferred scale: they can involve a few (commonly) or all (rarely) 
of the elements in the system. Their sizes and lifetimes follow power-law distributions, 
whose exponents are related. The numerical values of the exponents depend on the 
spatial dimensionality of the system, the spatial symmet ries of the l ocal forces and 
redistributive channe ls, the strength o f the local forces (iField et al.lll995l ). and the 
level of conservation ( Olami et al. 19921 ). % 



. Over the long term, the system tends to a critical state, which is stationary on average 
but not instantaneously. For example, on average, the power input by the external 
driver equals the energy per unit time released by avalanches. But there are fluctua- 
tions, because, at any instant, a random amount of energy is stored in metastable local 
reservoirs. 

Avalanche dynamics are generically observed in nature when conditio ns (i)-(iii) ar e 



met, and properties (iv)-(vi) emerge irrespective of the detailed microphysics (j Jensen! 1 19981 ). 
Likewise, in this paper, we remain agnostic about the microphysics of pulsar glitches; the 
statistical analysis presented below makes no assumptions in this regard. Nevertheless, it 
is striking that the traditional glitch paradigm — collective unpinning of quantized su- 
perfiuid vortices interacting with an inhomogeneous, slowly decelerating crust — conforms 



^In this respect, far-from-equilibrium critical systems differ from equilibrium critical systems (e.g. second- 
order phase transition in a ferromagnet), whose exponents depend only on the dimensionality of the system 
and its order parameter (s). 
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closely with (i)-(vi) (lAnderson &: Itohlll975l: lA. 
paradigm, based on crust fracture (jAlpar et al 



1996 



par et al.lll996l). So too does an alternative 



Middleditch et al.ll2006l ). whose terres- 



trial counterpart (pla te tectonics) is renowned as an archetype of self-organized criticality 
( ISornette et al.lll99ll ). We elucidate the analogy briefly before continuing. 



Consider a rectilinear array of quantized vortices, each carrying circulation k, spaced 
evenly according to Feynman's rule (47rz//K vortices per unit area) in the neutron superfiuid 
permeating the inner crust of a neutron star. A small percentage of the vortices are pinned to 
defects and/or nuclei at randorn loca tions in the crustal lattice, clustered to varying degrees 
(lAlpar et al.l Il996l : IWong et al.l I2OOII ) As v decreases gradually due to electromagnetic spin 
down, most vortices move apart and the outermost ones are expelled. However, the pinned 
vortices stay (nearly) fixed, in metastable reservoirs separated by relaxed zone s [see (ii)]. 



creeping slowly between adjacent pinning sites in response to thermal fluctuations (ILink et al. 



I993I ). The reservoirs are identical to t he capacitive eleme nts (v ortex traps s urrou nded by 
vortex depletion regions) postulated by lAlpar et al.l (119961 ) and I Wong et al.l (1200 ll ). They 
may be seeded by starquakes, which create large numbers of fresh lattice dislocations with 
deep pinning potentials, or they may emerge spontaneously in the self-organized critical state, 
as successive generations of vortex avalanches traverse the crust. As the pinned vortices 
increasingly lag the regular, unpinned array, a gradient in vortex density is established, and 
the local Magnus force on a pinned vortex rises. When the pinning threshold is overcome, 
a pinned vortex unpins and moves abruptly away from the pinning site (stick-slip motion), 
disturbing the local superfiuid velocity field (and hence the Magnus force) appreciably. Often, 
this is enough to push neighboring, barely subcritical, pinned vortices over their thresholds, 
triggering an avalanche. Most vortices in the avalanche rejoin the regular, unpinned array, 
and the crust spins up proportionately to compensate. The time-scale for a vortex to adjust 
locally to the Magnus and pinning forces is much shorter than z//z>, in keeping with (iii). 

Classic laboratory experiments on magnetic flux vortices in a type II superconductor 
(e.g. NbTi) immersed in a slowly c hanging magnetic field, an exactly analogous system, 
clearly exhibit properties (iv)-(vi) (IField et al.lll995l ). Vortices are expelled mostly in a 
continuous flow (cf. steady spin down) and occasionally in avalanches (cf. glitches). The 
distribution of avalanche sizes is measured to be a power law over several decades, whose 
exponent depends on the strength of the applied magnetic field (which controls the vortex 
spacing and hence the strength of the vortex- vortex interaction). The temporal fluctuation 
spectrum scales as an inverse power of frequency at high frequencies. After initial transients 
die away, the superconductor fluctuates around a self-organised critical state, called the Bean 
state, where the Lorentz force acting on each vortex is everywhere equal to the maximum 
pinning force. 
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If the pinning sites are sparsely distributed, so that global (mean-field) forces dominat e 
local forces between pinned vortex clusters, scale invariance breaks down (jjenseni Il998l ). 
Avalanches still occur, but they are distributed narrowly around a characteristic size and 
lifetime, involving nearly all the vortices instead of small, independent subsets. In this 
regime, avalanches recur quasiperiodically, not stochastically. Similar behavior is observed 
when the external driver acts too rapidly, but this situation never arises in pulsars. 

Scale invariant avala nche dynamics a nd self-organise d critical states are ob served widely 
elsewhere, in sandpiles (Bak et al. _ 1987 ). earthquakes ( Sornette et al.l E99ll ). so lar fiares 
( Lu fc Hamilton 1991 : Wheatland 200ol ) . and bursts from soft-gam ma-ray repeaters jGogiis et al, 
2000 ). The analogy with p ulsar glitches has been pointed out by iMorlev fc Garcfa-Pelavcj 



1993af) and ICarrollI (119981 ) and modeled using a cellular automaton by iMorley fc Schmid 



fll996h . 



Data 



Table [T] lists all 285 glitches discovered up t o the time of writing and known t o the 



authors. It is compiled from published sources ( Shemar &: Lvne 



Wan 



200 



et al.ll2000l:lHobbsll2002l : [Krawczvk et~aDl2003l : ljanssen fc Stappers 



1996 



Lvne et al 



2006 



2000 



Middleditch et al 



Peraltall2006f h the Australia Telescope National Facility Pulsar Catalogue ([Manchester et al.l 



20051 ). which can be accessed on-line at http: //www. atnf . csiro . au/research/pulsar/psrcat 



and unpublished data communicated privately by M. Kramer, D. Lewis, and A. G. Lyne. For 
each pulsar, the table lists its J2000 coordinates, and the number of glitches detected (Ng). 
The earliest and latest epochs observed (tmm and tmax respectively) are recorded separately in 
Table [2] for the nine pulsars with iVg > 5. An asterisk signifies that segmented data spans are 
not specified in the cited references; in this situation, tmin and tmax are estimated by eye from 
spin-down histories graphed in the cited references, where available, or else from the first 
and last glitches by default. For each glitch. Table [1] lists its epoch, the fractional increase 
in spin frequency Au/u, and one or more bibliographic references. Uncertainties are quoted 
as a trailing integer in parentheses, corresponding to an absolute number of days for t [e.g. 
MJD 51141(248) means MJD 51141 ±248] and an uncertainty in the last signifcant digit for 
Au/u [e.g. 0.04(2) means 0.04 ± 0.02]. For some newly discovered glitches, the information 
is incomplete. Epochs and sizes have been measured for 271 and 250 glitches respectively. 
Other parameters, like the healing fraction and post-gl itch relaxation time-scale, are omitted 
as they are not analyzed in this paper; please consult iPeraltal (120061 ) and references therein 
for a full catalog. 
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4. Size distribution 

4.1. Scale invariance 

If pulsar glitches are the result of an avalanche process, their size distribution should be 
scale invariant in any individual pulsar, with probability density function 

p{Au/u) oc {Au/u)~^ . (1) 

The exponent a is set by the dimensionality [fl and symmetries of the local forces, which are 
likely to be universal, and the strength and level of conservation of these forces, which are 
functions of temperature and therefore not universal (see One therefore expects a to 
differ from pulsar to pulsar. As a corollary, the aggregate size distribution drawn from all 
pulsars is not expected to be a simple power law of the form ([1]). 

To test these ideas, we construct the observed cumulative size distributions of the nine 
known pulsars with Ng > 5. The selection criterion A^^g > 5 is arbitrary; it seeks to limit the 
impact of random errors while testing as many objects from Table [1] as possible. We then 
compare the data against the theoretical cumulative distribution 

_ (A./.)-'-(A./.).^; 



/mm 



derived from ([T]). The theoretical distribution is normalized after restricting it to the domain 
(Az//z/)min < Au/u < (Az//z/)max, wlicre (Az//z/)min and (Az//z/)max are the smallest and 
largest glitches observed in that pulsar respectively, quoted in Table [H There are more 
sophisticated ways to choose (Az//z/)min and (Az//z/)inax7 which we consider further below, 
but this is a conservative starting point. 

For each object, we choose a to minimize the Kolmogorov-Smirnov (K-S) statistic D, i.e. 
the maximum unsigned distance between the curves. The numerical results are recorded in 
Table [3|, while the measured and theoretical cumulative distributions are plotted together in 
Figured! (Cumulative distributions are free of binning bias.) The goodness of the fit at the 
optimal value of a is characterized by Pks, defined such that 1 — Pks equals the probability 
that the K-S null hypothesis (that the two data sets are drawn from the same underlying 
distribution) is false. The 1-cr lower and upper bounds a_ and a+ mark the range of a 



■^The effective dimension need not equal three. For example, it may equal two in a rectilinear vortex array 
or faulting crust, where t he local forces act in the transverse plane, and three in a turbulent vortex tangle 
|Peralta et al.ll2005l . l20oi ). 

*The K-S test is most sensitive to discrepancies near the median. An alternative test, based on Kuiper's 
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where the null hypothesis is rejected with less than 68% confidence. Note that the interval 
[a_, a^] is asymmetric about the optimal a and widest for the best fits. 

The results in Table [3] confirm what is apparent by eye from Figured! the null hypothesis 
that the size distribution is described by a power law for all nine pulsars with A^g > 5 is 
not ruled out at the 1-a level of confidence. In turn, this is consistent with the avalanche 
hypothesis. However, in two objects, namely PSR J0537— 6910 and PSR J0835— 4510, the 
agreement is marginal. Interestingly, these two objects are also the only ones discovered so 
far that are believed to glitch quasiperiodically (iLyne et al.lll996l : iMiddleditch et al.ll2006l ). 



4.2. PSR J0537-6910 and PSR J0835-4510 

Quasiperiodicity is a natural feature of avalanche dynamics when mean-field forces over- 
whelm local interactions, as described in §21 We explore its manifestation in glitch waiting 
times in §3 With respect to glitch sizes, we note that avalanche s in the quasi periodic regime 



tend to be distributed narrowly around a characteristic size (jJensenl 119981 ). This can be 



modeled crudely by adding a term proportional to (5[Ai//z/ — (Az//z/)c] to ([T]), viz. 

p{Au/u) = a{Au/u)-'' + CA^u/u - (Au/u),] , (3) 

where (Az//i/)c denotes the characteristic size, and the scale invariant and quasiperiodic 
components are weighted by the constants Cs and Cp respectively. Normalization fixes Cg 
in terms of Cp (or vice versa), with Cg + Cp 7^ 1 in general. The associated cumulative size 
distribution is given by 

P{Av/v) = CpH[Av/v - {Av/v)^] 

{l-Cp)[{Au/uf-'^~{Au/u)]^^] 



' mm 



where H{-) denotes the Heaviside step function. 

Parameters determined by fitting (jl]) to the data are recorded in Table HI while the 
corresponding measured and theoretical cumulative distributions are plotted together in 
Figure [2l The fits are much improved, with Cp ~ 0.2 in both objects — although, to 



statistic, mitigates this bias (|Press et al.lll986f ). It will be worth implementing when more data become 
available. The K-S test is also inefficient if the underlying probability density function contains a narrow 
notch, where the probability vanishes. Again, there is insufficient data at hand to look for such a notch; it is 
difficult to find in a cumulative distribution, and the probability density function is biased by binning when 
N„ is small. 
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be fair, the delta-distributed component is not strictly required, at least not at the 1-a 
level. Importantly, the delta-distributed component contains only ~ 20% of the glitches, 
not all of them. This is consistent w ith the historical interpretation of the pulsar data 
(ILyne et al.lll995l : [Marshall et al.ll2004j ). It is also seen in self-organized critical systems like 
sandpiles, where large, system-spanning, quasiperiodic avalanches of a characteristic size 
are interspersed with small, random ly timed avalanches, which are power-law distributed 
jRosendahl et al.lll993l : Ijensenlll998[ l. 



4.3. Upper and lower cut-offs 

Strictly speaking, it is incorrect to normalize P^Au/u) by choosing (Ai//i/)jnm and 
(Az//z/)max to bc the smallest and largest glitches observed in a pulsar respectively. A better 
choice of (Atz/z/l min is the actual resolution of the timing experiment, which varies with 



object and epoch. iJanssen fc Stapperd (120061 ) simulated detection of a microglitch in a noisy 
time series and obtained {Ai>/i>)^[n = 1 x 10~^^. Usually, this information is not provided 
explicitly and must be estimated from the size uncertainties quoted for detected glitches. 
All the same, the smallest glitch observed is likely to be a reasonable estimate of Az/min, 
because the occurrence probability increases steeply as Az/ decreases, according to Table [31 
On the other hand, (Az//z/)jnax is limited by the total observing time. Its true value exceeds 
the largest glitch observed, but not by much, because the occurrence probability decreases 
steeply as Ai/ increases. 

To quantify these effects, we allow (Az//z/)inin to vary between 0.5 and 1.0 times the 
smallest glitch size observed, (Az//i/)inax to vary between 1.0 and 2.0 times the largest glitch 
size observed, and fit equation (jl]) again to the data. For every object, (Az//z/)min and 
(Az//z/)max shift ouly slightly, and a stays within the range [a_,a+] in Tables [3] and HI This 
confirms that the smallest and largest glitches provide reasonable estimates of (Az//i/)inin 
and (Az//z/)max- At the 1-a level, the constrained and unconstrained fits are both consistent 
with the data. 



4.4. Aggregate distribution 

Figure [3] displays the cumulative size distribution for the glitch population in aggre- 
gate, together with the best power-law fit of the form ([2]). The fit is poor. When all 250 
glitches with measured sizes are included, the best fit corresponds to a = 0.96, (Az//z/)min = 
9.5 X 10-1^ (Az//z/)rnax = 2.0 X 10-^ and Pks = 7.1 x 10"! When the glitches from 
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PSR J0537-6910 and PSR J0835-4510 are excluded, the best fit corresponds to a = 0.98 
and Pks = 3.2 x 10~^, with (Az//z/)jnm and (Az//z/)max as before. Either way, we can state 
confidently that the theoretical and observed data are drawn from different underlying distri- 
butions. This is not surprising; the results in Figure [T] and Table [3] demonstrate clearly that 
the size distribution in individual pulsars is consistent with being scale invariant, but that a 
differs from object to object. Hence the aggregate distribution is expected to be a weighted 
sum of power laws, not a pure power law. Accordingly, the size distribution i n individual 



pulsa rs is a more direct probe of glitch physics than the aggregate distribution (ILyne et al. 



20001 ) . The aggregate distribution can be inverted, in principle, to determine how a is dis- 
tributed across the pulsar population. We defer this exercise until better historical estimates 
of (Az//i/)inin and (Ai//i/)inax, and more data, become available. 



Janssen fc Stapperd (120061 ) claimed that the glitches in PSR J1740— 3015 are drawn 
from a flat size distribution in log(Ai//z/), i.e. a = 1, with Pks = 0.902. This agrees with 
the results in Table [3l 



Lyne et al.l (120001 ) noted some evidence for an excess of large glitches, which is cor- 
roborated to some extent by Figure [3^. However, the excess largely disappears when the 
quasiperiodic glitchers are excluded, as in Figure [3]d. Large glitches do not originate prefer- 
entially from any particular class of object. While it is true that the most active objects (e.g. 
PSR J0537— 6910 and PSR J0835— 4510) experience relatively large and narrowly distributed 
glitches, with Az//z/ > 10"'', other active objects (e.g. PSR J1740— 3015) experience a mix 
of small and large events, and there are several objects (e.g. PSR J1806— 4212) which have 
only glitched once, with Av/v > 10~^ for that single glitch. Furthermore, although PSR 
J0534+2200 is sometimes portrayed as unusual for not experiencing large glitches, its size 
distribution is actually relatively flat (a ~ 1.2). There is every reason to expect that it will 
experience large glitches in the future, but these events will be slow in coming, because PSR 
J0534+2200 builds up differential rotation between the crust and superfiuid at a relatively 
slow rate, as we show in ^ § 



^ IWong et al.l (j200ll ) argued that the relative angular acceleration of the crust and superfiuid in the Crab, 
inferred from the activity parameter, is much smaller than expected given the large Ajy/v ~ 10^"* observed 
during glitches. This paradox is resolved if most of the differential rotation is being stored temporarily, in 
advance of a large glitch in the future. 
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5. Waiting-time distribution 



5.1. Poisson process 

If pulsar glitches are the result of an avalanche process, they should be statistically 
independent events. To understand why, recall that a system in a self-organized critical 
state configures itself into many metastable stress reservoirs insulated by relaxed zones ([|2]). 
Every avalanche relaxes one reservoir, typically occupying a small fraction of the system, and 
the next avalanche occurs at random, typically far from its predecessor. There is essentially 
no interfer ence between successive avalanches; this is verified empirically in tests with cellular 
automata (jJenserull998l ). Avalanches in the tail of the size distribution, which relax the whole 
system, are an (extremely rare) exception. 

Given statistical independence, and assuming that the system is driven at a constant 
(mean) rate, the avalanche model predicts that the tirae between successive g litches, At, 
termed the waiting time, obeys Poisson statistics (jjensen 1998 : Wheatland 2000 ). Hence, in 
any individual pulsar, the waiting-time distribution is exponential, with probability density 
function 

p(A, At) = Aexp(-AAt) . (5) 

The mean glitching rate A is different for every pulsar. It depends on the rate at which 
differential rotation builds up between the superfiuid and the crist (oc z>) as well as the 
capacity to store the differential rotation (e.g. strength of pinning, rate of vortex creep, shear 
modulus of the crust). The storage capacity is presumably controlled by thermodynamic 
variables like temperature, as well as the inhomogeneous nuclear structure of the crust. We 
do not expect A to change appreciably during four decades of pulsar timing. In principle, 
however, as more data are collected in future, this claim can be tested by using a Bayesian 
blocks algorithr n to divide th e time series into a sequence of Poisson processes Av ith piecewise- 
constant rates ( Scargle 1998 : Wheatland 2000 : Connors fc Carraminana 20031 ). 



The avalanche model makes a further powerful prediction. Suppose the system tends 
to a stationary, self-organised critical state, in which global driving balances local release in 
a time-averaged sense [property (vi), §2], i.e. there is no secular accumulation or leakage of 
stress. Stationarity implies that the mean waiting time (At) = A~^, multiplied by the rate at 
which crust-superfluid differential rotation builds up (ez>), equals the mean ghtch size (Au), 



^ It is tempting to assume that A is constant over decades, because the thermodynamic variables that 
control storage capacity (e.g. temperature) are nearly constant on such a time-scale. Yet the Sun provides a 
cautionary counterexample: the dynamics of subp hotospheric turbulence, and hence the rate of solar flaring, 
vary with the 11-yr solar cycle ( WheatlandluOOOl ). 
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I.e. 



A = -ez>/(Az/) 



(6) 



Here, 27rez> is the relative angular acceleration of the crust and superflu id. Importantly, 



glitch data allow e to be measured directly in principle (jWong et al. 



200 ll ). However, there 



is a serious question as to whether stationarity is achieved in practice, during the 40 yr that 
a typical pulsar has been observed. We discuss this issue further below. 

To test the above ideas, we compare the observed cumulative waiting-time distributions 
of the same nine pulsars as in §11 with A^^g > 5, against the theoretical cumulative distribution. 
In order to make the comparison fairly, we must first adjust for observational selection 
effects. Any given obesrvation can detect waiting times in a range Atmin < At < Atmax- 
The upper limit Atmax is set by the total data span available for that pulsar, i.e. Atmax = 
^max — ^min- The lower limit Atmin is different at different epochs. It is set by the gap 
between data spans in which a glitch is localized. For small glitches, the glitch epoch is 
determined by requiring continuity of pulse phase across the glitch. For larger glitches, 
where the phase winding n umber is arabiguo us, the epoch is taken to be halfway between 
the bounding observations (I Wang et al.ll2000l ). Either way, At^[^ is differe nt for each glitch , 



and is twice the absolute value of the epocha l uncertainty quoted in T able [T] (ILyne et al.ll2000 



Wang et al.ll2000l : I.Tanssen fc Stappersll2006l : iMiddleditch et al.ll200fil ). Let /(Atmm)c^(Atmin) 
be the observing-time- weighted probability that, when a glitch occurs, Atmin hes in the 
range [At^^i^, At^[^ + d{At^[^)], and let the smallest and largest values of At^[^ be At^^]^ 
and At[j^^ respectively, recorded in Table \5\ for the nine pulsars with A^^g > 5. Then the 
cumulative waiting-time distribution is given by 



P(A,At) 



At 



(» 



At 



«) 

min 
Atn 



d(ACJ/(ACJ 



X / d{At') p{\, At') 



(7) 



1 



E 



exp(— AAtjnin) — exp(— AAt) 



At„ 



=At 



«) 



exp(— AAt„ 



exp(-AAi(:max) 



where each glitch is weighted equally in the sum in ([8]) as a first approximation. 

In Figure HJ we plot as cumulative histograms the measured waiting-time distributions 
of the nine pulsars in Figure [TJ The theoretical curves ([8]) are overlaid, with Atmin and 
Atmax chosen according to the second through fourth columns in Table |5l For each object, 
we choose A to minimize the K-S statistic D. The fitting parameters are displayed in the 
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fifth through eighth columns in Table [51 As in §11 the goodness of the fit is characterized by 
the K-S probability Prs? with Pks < 0.32 in the interval [A_, A_|_]. 

For all nine pulsars in Figure H] and Table 0, the null hypothesis that the waiting-time 
distribution is described by Poisson statistics is not ruled out at the 1-a level. The data are 
therefore consisten t with an ava l anche process. An exponential waiting-time distribution was 



first postulated by lWong et al.l ( l200ll ) for PSR J0534+2200, based on timing data up to and 
including the glitch on MJD 51452. These authors obtained A = 0.53 yr~^ and Pks = 0.7, 
marginally outside the 1-a range in Table [5l The data analyzed here confirm that waiting 
times are consistent with Poisson statistics in several glitching pulsars, affording a key insight 
into the physics of the glitch mechanism. The implications of this result are discussed further 
in m 



5.2. Quasiperiodicity 

For seven pulsars in Figure HI the Poisson distribution affords an excellent fit, both 
formally and by eye. However, for PSR J0537— 6910 and PSR J0835— 4510, the fits are 



margi nal at the 1-a level. [Indeed, in an earlier analysis of PSR J0835— 4510, IWong et al. 



(120011 ) excluded a Poisson distribution with 96% confidence, on the basis of fewer data.] 
These are the same two objects whose size distributions are exceptional, and which are 
observed to glitch quasiperiodically. 

Taking the same approach as in ^ we model the quasiperiodicity crudely by adding a 
periodic component to viz. 

p{X, At) = C'^X exp(-AAt) + C'p6{At - At,) . (9) 

In ([9]), Cg and C'^ are the normalized relative weights of the Poisson and periodic components 
respectively, and At, is the characteristic period. The associated cumulative distribution, 
weighted by At^in, is obtained by substituting ([9]) into ([7]), yielding 



1 r 

P(A,At) = ^ 5^ lc'^H{At-At,) 

^ At ■ -At«' 

(1 - Cl)[exp{-XAtmin) - exp(-AAt)] 



exp(— AAtmin) — exp(— AAt 



max ) 



(10) 



The two-component model ([9]) yields improved fits to the data, with ~ 0.25 for 
both objects. The fits are graphed with the data in Figure [5|, and the best-fit parameters 
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are recorded in Table El We obtain At^ = (0.3 ± 0.1) yr and Ate = (2.8 ± 0.1) yr for PSR 



J053 7 — 6910 and PSR J0835— 4510 respectively, in accord with previous authors (ILyne et al. 
19961 : iMiddleditch et al.l 120061 ). Significantly, the data imply Cp ~ C' In other words 



the delta-distributed component accounts for the same fraction of the size and waiting-time 
distributions, even though the sizes and waiting times are statistically independent. This 
raises confidence in the model and suggests that a quasiperiodic component is indeed present 
and distinct. It also suggests that the quasiperiodic component coexists with the Poisson 
component, instead of completely displacing it. Vela, for example, is likely to possess an ex- 
tensively connected network of capacitive element s, but it also conta ins smaller subnetworks 
that are disconnected from the mai n network; cf. lAlpar et al.l (19961). T his is natural for an 



avalanche process, as noted in §4.21 (IRosendahl et al.l 



1993 



Jensenlll998r ) 



5.3. Mean rate 

Stationarity of the avalanche model over long time intervals implies a relation between 
the Poisson rate, driving rate, and mean glitch size, given by ([6]). Unfortunately, for a < 2, 
(Az/) is dominated by large glitches near the upper cut-off of P(Az//z/): 



(Az.) 



a - 1 



a-2 



Az^upper 

/Azv 
Az/, 



upper 



a-1 



Az/, 



upper 



upper 



a > 2 
1 < a < 2 
a < 1. 



In ( ITT]) . Az/iower and Az/^pper are the physical lower and upper cut-offs of the probability 
distribution function ([T]). As large glitches are rare, stationarity is not achieved during the 
40-yr interval over which a typical pulsar is observed; the largest observed size, (Az//z/)max, 
cannot be equated reliably with the maximum size allowed physically. Likewise, (Az/) is 
approximated poorly by the mean of the observed glitches. In practice, therefore, it is 
impossible to estimate (Az/) without much longer monitoring. 

Physically, ez/ is the rate at which differential rotation builds up between the crust 
and superfiuid. Hence, in the vortex unpinning model, e gives the time-averaged fraction of 
pinned vortices or capacitive elements. We can use equation (fTT!) to place limits on e, at least 
in principle. For example, the inequalities Az/iower < Az/min and Az/max < Az/upper < must 



^For 1 < a < 2, (Av) is dominated by the upper cut-off, but the normalization of p(Az//z/) is dominated 
by the lower cut-off. 

®The error in Ai/ is d{Ai') /da multiplied by the error in a. 
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always be satisfied. For PSR J0358+5413, assuming a = 2.4, we find (Az/)/z/ < 1.1 x 10"^° 
and hence e < 7 x 10~^. This is lower than previous estimates of the pinned fraction for 
objects of that age, but in line with previous estimates of the pinned fraction in young 
objects like the Crab (ILyne et al.ll2000l : IWong et al.ll200ll ). For the remaining eight objects, 
the above inequalities lead to upper limits on e which are greater than unity and hence 
not useful. As a crude experiment, we check the result of setting Ai/upper = 2 x 10~^, the 
largest glitch observed in any pulsar over the last 40 years, for every pulsar. We obtain five 
slightly more useful upper limits (e < 4 x IQ-^, 0.2, 0.8, 0.8, and 0.7 for PSR J0534+2200, 
PSR J0631+1036, PSR J0835-4510, PSR J1341-6220, and PSR J1740-3015 respectively). 
However, we emphasize that these values are still problematic, because there is no guarantee 
that a total effective observation interval of 40 yr x 101 pulsars is long enough in aggregate 
for stationarity to be observed. Moreover, these values are based on assuming that all pulsars 
have the same physical Ai/upper? which is not necessarily the case. 

Figure [6] displays the cumulative histogram of Poisson rates derived from the best-fit 
waiting-time distributions in Figures H] and [51 Let g(A) denote the rate probability density 
function, such that q{X)dX is the probability that the mean rate lies in the interval [A, X + dX] 
in a given pulsar. There is no obvious theoretical reason to prefer a particular analytic form of 
g(A), which is controlled by the physics of the global driver, not the scale invariant avalanche 
dynamics. In addition, the data in Figure [6] are insufficient to specify the analytic form of 
g(A) uniquely. However, motivated by the rate distribut ion observed in solar fiares, which 
is measured reliably to be exponential (IWheatlandl l2000l ) . we find that a distribution of the 
form 

g(A) = (A)-iexp(-A/(A)) (12) 



l.StolyT ^ including quasiperiodic glitchers {left 



fits the data satisfactorily, with (A) 
panel) or (A) = 1.2l^Q'4yr~^ excluding quasiperiodic glitchers (right panel). Formally, the 
K-S probabilities are Pks = 0.9946 and 0.82 respectively. 

The distribution is incompletely sampled below an effective minimum rate Amin, which 
is set by Atmax- To illustrate, if we proclaim five glitches (arbitrarily) to be the minimum 
number needed for a reliable determination of A, we obtain Amin = 5At~?^^ ~ 0.2 yr^^. 
Careful modeling of this observational bias is deferred to a future paper. We describe a first 
attempt in §5.41 



5.4. Aggregate distribution 



The nine pulsars in Figure H] account for only 108 out of a total of 285 observed glitches. 
Most glitching pulsars have only glitched once or twice, but they still contribute statistical 
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information regarding waiting times, the former via lower limits on At. While these data 
cannot usefully constrain P(A, At) in individual pulsars, they feed into the aggregate waiting- 
time distribution and hence constrain q{X) more tightly than in §5.3[ 

In Figure [TJ we present the aggregate waiting-time distribution Pagg(At) including {left 
panel) and excluding [right panel) the quasiperiodic glitchers PSR J0537— 6910 and PSR 
J0835— 4510. The histogram is constructed to include all 182 waiting times in those objects 
that have glitched more than once. The K-S test confirms that the aggregate distribution is 
fitted poorly by a single, constant-rate, Poisson distribution of the form Furthermore, 
when we weight ([5]) by the exponential rate distribution f|T2|) . |^ viz. 

Ai POO 

Pagg(At) = / d{At') / dX'qiX')piX',At'), (13) 
Jo Jo 

the fit remains poor. For example, the dotted curves in Figure [7] are computed by evaluating 
( IT3l) with the mean values (A) = 1.3logyr^^ (left panel) and (A) = 1.2lo;4yr~"'^ (right panel) 
extracted from Figure El They yield Pks = 4.3 x 10^^ and 2.4 x 10^^ respectively. If, 
instead, we adjust (A) to maximize Pks while fitting f|T2l) and f|T3|) to the observed Pagg(At), 
as shown by the dashed curves in Figure [71 we obtain (A) = 1.1 yr~^, Pks = 0.18 (left panel) 
and (A) = 0.92 yr~^, Pks = 0.31 (right panel) respectively. 

We can exploit the extra information in Figure [71 from objects with 2 < Ai'g < 5 to 
determine q(X) more accurately. To do this, we assume a rate probability density function 
of the form 

gA = {Xy'H(X - A^in) exp[-(A - A^i„)/(A)], (14) 
normalize P(A, At) over the range [Atmin, Atmax], and evaluate f[T^ to obtain 



POO 

Pagg(At) = (A)-i / dX' exp[-(A' 



A„i.)/{A)1 



X (15) 

l-exp(-A'At, ' ^ ^ 



^max / 



In ([T5|l . we neglect for simplicity the observational bias introduced by uncertainties in glitch 
epoch discussed in §5. It lacking fuller information, we take Atmin = and Atmax = 28.03 yr 
for all pulsars. Excellent fits are obtained using f[T5l) . We find (A) = 0.54 yr~\ Amin = 0.25 
yr~^, and Pks = 0.76 for all nine pulsars with A^g > 5, and (A) = 0.43 yr^^, Amin = 0.25 yr~^. 



^We compute Pagg(Ai) theoretically as a weighted sum of independent Poisson processes. In the same 
way, the waiting-time distribution for decays observed from a mixture of radioactive isotopes is a sum of 
constant-rate Poisson distributions, one per isotope, weighted by isotopic abundance. 
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and Pks = 0.98 when the quasiperiodic ghtchers are excluded. The fits are plotted as solid 
curves in Figure [7] ( /e/t and right panels respectively). We verify the results by referring back 
to the measured g(A) distribution. Substituting the fitted values of (A) and Amin into ffT^ . 
we obtain the solid curves in Figure El with Pks = 0.52 {left panel) and 0.25 {right panel) 
respectively. Alternatively, the previous fits (A) = 1.1 yr~^, Amin = and (A) = 0.54 yr~^, 
Amin = yield Pks = 0.72 {left panel, dashed curve) and 0.52 {right panel, dashed curve). 

In all cases, Pagg(At) points to the existence of more low-rate objects than the A'g > 5 
sample in Figure E] predicts. Specifically, up to ~ 30 % of the population of glitching pulsars 
can have A < 0.25 yr~^ while still reproducing Pagg(At). We emphasize again that Figure [7] 
constrains q{X) more tightly than Figure O because it contains information about At from 
1.3 times as many glitches, including useful information from pulsars which have glitched 
twice. 

Undetectable microghtches probably occur between detected glitches without our knowl- 
edge, given that p(Az//z/) is scale invariant. This effect subtracts from the lower end of the 
At distribution and adds to the upper end. We do not correct for it here, because it is hard 
to quantify without better statistics. On two occasions, a pair of g litches occurred on th e 
same date, once in the same pulsar, and once in different pulsars faramer fc Lyne 2005). 



We take At = for these pairs. Phase connected timing mitigates duty cycle biases, but it 
does not eliminate them. 



5.5. Fluctuation spectrum 

The temporal fluctuations in a stochastic signal x{t) carry independent statistical in- 
formation about the underlying physical process. The power spectrum S{f), where / de- 
notes the Fourier frequency, is related to the temporal autocorrelation function G{t) = 
{x{t)x{t + T)) — {x{t))'^ (where the average (...) is performed over t for a stationary process) 
through the cosine transform 

POO 

S{f) = 2 dTG{T)cos{2nfr) . (16) 
Jo 

In general, for an avalanche process, the power spect rum depends jointly on the size. 



waiting-time, and lifetime distributions of the avalanches (lJensenlll998l ). For glitches, how- 
ever, the lifetimes are too short to measure with current technology (see ^HD- If, furthermore, 
we restrict attention to the unit- impulse signal x{t) = J2i ^{t~U), where tj denotes the epoch 
of the i-th glitch, then the sizes drop out of the problem too. The power spectrum then car- 
ries exactly the same information as the waiting-time distributions P(A, At) and Pagg(At), 
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with 

for any individual pulsar, and 

for the pulsar population in aggregate. 

At high frequencies / ^ (A), equations (ITTl) and (ITSl) [with g(A) given by (fT2l) ] scale 
as /~^, with 0{f~^) and 0(/~^sin2/) corrections. These scalings are modified if the delta 
function in x{t) is replaced by a nonsingular window function that embodies the shape of 
the signal from an individual avalanche. It will be instructive to revisit this question when it 
becomes possible to resolve the lifetimes of individual avalanches, e.g. in single- or giant-pulse 
timing experiments. 



6. Discussion 

In this paper, we analyze the size and waiting-time distributions of pulsar glitches, 
taking advantage of the enlarged data set produced by the Parkes Multibeam Survey. We 
conclude that the data are consistent with the hypothesis that pulsar glitches arise from an 
avalanche process. In each of seven pulsars with A'^g > 5, the size distribution is consistent 
with being scale invariant across the observed range of Az/ (up to four decades), and the 
waiting-time distribution is consistent with being Poissonian. These features are natural if 
the system is driven globally at a constant rate (as the pulsar spins down), and each glitch 
corresponds to a locally collective, threshold activated relaxation of one of the many spatially 
independent, metastable stress reservoirs in the system (e.g. via a vortex unpinning or crust 
cracking avalanche). In two pulsars, PSR J0537— 6910 and PSR J0835— 4510, the dynamics 
may include a second, quasiperiodic component, comprising ~ 20% of the events. The 
size and waiting-time distributions of the quasiperiodic component are narrowly peaked, as 
expected for rare, system-spanning avalanches, which relax a large fraction of the total stress 
accumulated in the system. This two-component behavior is observed widely in self-organised 
critical systems, including experiments on magnetic flux vortices in type I I superconductors. 



which are closely analogous to neutron star superfluids (IField et al.lll995l ) 



The power-law exponent of the size probability density function differs from pulsar to 
pulsar, spanning the range —0.13 < a < 2.4. Calculating a theoretically from first principles 
is a deep problem which has not yet been solved for other self-organised critical systems. 
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let alone glitching pulsars, although some prog ress has been made on two-dimensional sand- 
piles using renormalization group techniques (jPietronero et al.lll994j : |Jensenlll998l ). In the 
mean-field approximation, which is exact in four or more dimensions, theoretical calcula- 
tions on sandpiles (and other systems in their uni versality class ) yield a = 1.5, whereas 
three-dimensional cellular automata output a = 1.3 (jJensenlll998l ). 



The size distribution transmits two important lessons concerning the microphysics of 
glitches. First, the fact that a differs from pulsar to pulsar implies that the strength and level 



of conservation o f the local (e.g. pinning and intervortex) forces also differs (jOlami et al.lll992 



Field et al.lll995l ). By contrast, in equilibrium critical systems like ferro magnets, a depends 
only on the dimensionality of the system and its order parameter and is therefore universal. 
Second, except for the two pulsars which show quasiperiodicity, a appears to vary smoothly 
with spin-down age, with a ~ 1.2 for the youngest pulsars (e.g. the Crab). Figure [3 depicts 
the trend between a and Tc. It is suggestive; after all, local pinning forces do depend on 
temperature and hence Tc. Interestingly, however, there is no clear trend between a and u, 
even though the mean vortex spacing (and hence intervortex force) is proportional to z/^/^. 
It will pay to study these trends more thoroughly as more glitch data is collected. 

An avalanche process predicts a specific relation between the distributions of glitch sizes 
Au and lifetimes T (as opposed to waiting times At). Specifically, in a self-organized critical 
state, the lifetime probability density function is also a power law, p{T) oc T~^, with 



b=l + {a - 1)72/73 



(19) 



The constants 72 and 73 are defined such that the cardinality of an a valanche scal es with 
its linear extent (L) as L'^^ and its lifetime (i.e. duration) scales as L'^^ (|Jensenlll998l ). Both 
72 and 73 depend on the effective dimensionality of the local forces and can be calculated 
numerically using a cellular automaton. In two dimensions, avalanches are compact, not 
fractal, and one has 72 = 2; in three dimensions, one has 2 < 72 < 3. At present, radio timing 
experiments cannot measu re T; most glitches are d etected as unresolved, discontinuous, spin- 
up events with T < 120s ( McCulloch et al. 1990 ). In the future, however, single- and/or 
giant-pulse timing experiments with more sensitive instruments (e.g. the Square Kilometer 
Array) will test this prediction. If confirmed, it will independently corroborate the avalanche 
hypothesis. 



^°In the Crab, some spin-up events seem to be resolv ed, e.g. at epochs MJD 50260 (T ~ 0.5 d) and MJD 
50489 (secondary spin up, T « 2d) (jWong et al.l 120011 ). If these are rare but otherwise standard ghtches 
originating from the long-T tail of the lifetime distribution, it is puzzling that other, shorter, but still resolved 
(and presumably more common) spin-up events are not observed, e.g. with T ^ 0.1 d or 0.01 d. Alternatively, 
the events at MJD 50620 and MJD 50489 may have been triggered by a different physical mechanism. 
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The mean glitching rates of the nine pulsars studied here are fairly narrowly distributed, 
spanning the range 0.35 yr^^ < A < 2.6 yr~^. The probability d ensity function fo r A is 
adequately fitted by an exponential, as for solar flare avalanches ( IWheatlandl l2000l ) . with 
(A) = l-S^Qg yr~^, or by an exponential with a lower cutoff, at Amin ~ 0.25 yr~^. A 
theoretical derivation of (A) from first principles is currently lacking, although estimates of 
how long it takes to crack the crust locally predict r easonable rates, if the critical strain 



angle approaches that of imperfect terrestrial metals (lAlpar et al.l 1 19961 : iMiddleditch et al. 

mm- 



Figure [TO] plots A versus Tc for the nine pulsars examined individually in this paper. 
There is no significant trend. The d ata are consistent with the notion that old pulsars glitch 
less frequently than young pulsars (IShemar &: Lynelll996l ). but they are equally consistent 
with the notion that the glitching rate is independent of age. 

Many authors have searched for a correlation between waiting time and the size of 
the next gli t ch. S uch a correlation appears to b e abse nt from the data, e.g. Figure 17 in 
Wang et al.l (120001 ) and Figure 10 in IWong et al.l (120011 ). At first blush, this is surprising: 
the vortex unpinning and crust fracture paradigms, which are driven by the accumulation 
of differential rotation and mechanical stress respectively, seem to be natural candidates for 
a 'reservoir effect'. Avalanche dynamics resolves this apparent paradox. The reservoir effect 
does operate locally, but the star contains many reservoirs, insulated from each other by 
relaxed zones, whose storage capacities evolve stochastically in response to the slow driver 
and avalanche history. During a glitch, a single reservoir (often small but sometimes large) 
relaxes at random via an avalanche, releasing its stored Au (and destabilizing neighboring 
reservoirs in preparation for the next glitch). Some of the Au is accumulated since the pre- 
vious glitch, but the remainder is 'borrowed' from earlier epochs, when some other reservoir 
relaxed instead. All self-organized critical syste ms share thes e dynamics; the waiting time is 
uncorrelated with the size of the next avalanche (jJenserul 19981 ). The only exceptions are large, 
system-spanning avalanches, which always have roughly the same sizes and waiting times, 
and which account for ~ 20% of the glitches in PSR J0537-6910 and PSR J0835-4510. 

A corollary of the previous paragraph is that the total Au released in glitches up to 
some epoch is less than the total crust-superfluid differential rotation accumulated since 
that epoch, viz. 

^Ai/i <e|z>|^Ati, (20) 
1=1 1=1 

where ez> is the relative angular acceleration of the crust and superfluid due to electromagnetic 
spin down. The 'stair case' described by (l20l) has been noted previously (jShemar fc Lyne 
19961 : iLyne et al.ll2000l ). both in quasiperiodic glitchers like PSR J0537— 6910 [e.g. Figure 8 
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in iMiddleditch et al.l (120061 )] . where the reservoir eff e ct is obvious, and in Poisson ghtchers 
hke PSR J0534+2200, [e.g. Figure 12 in IWong et al.l (|200ll )]. where the trend is more subtle 
because it reverts to the mean over long times, not after every glitch. Upon dividing fl2Ul) by 
Ng, and averaging over long times, the inequality becomes an equality (provided there is no 
secular accumulation of differential rotation in the system) and we recover ([6]). 

It is fundamentally impossible to measure e for individual pulsars with current data, 
because (Az/) is dominated by large (and therefore rare) glitches for a < 2. It is there- 
fore wrong to assume stationarity over a typical, 40-yr observation interval. Consequently, 
we are prompted to r eassess the familiar correlation between acti yity and spi r i-dow n age 



( IShemar fc Lynelll996l ). Our definition of ez> is id entical to i^M r- i i in iLvne et al.l (120001 ) (but 



for individual objects, not in aggregate) a nd Ap. inlWong et al.l (120011) . It is closely related to 
the original activity parameter defined by McKenna &: Lyne~(jl990 ). which equals NgU^^eli/]. 
For PSR J03 58+5413, we meas ure e < 7 x 10^^, lower than the aggregate value 0.017± 0.002 
measured by iLyne et al.l (120001 ) for objects with Tc > lOkyr (binned by semi-decades in z>). 
Interpreted in terms of the vortex unpinning model, this result suggests that 0.007-2 % 
of the angular momentum outfiow during spin down may be stored in metastable reservoirs 
on average over time. On the other hand, five other objects have 0.04 < e^ax < 0.8, under 
the questionable assumption that the maximum physical size is Az/, 



upper 



2 X 10"^ in all 



pulsars. Our data are therefore inadequate to update usefully the value Ag/|z^| 
measured by Iwong et all J200lh for PSR J0534+2200. 



1 X 10" 



In the context o f vortex unpinning, it has been argued that the aggregated e measured 
by iLyne et al.l (120001 ) partly corroborates the hypothesis that younger pulsars are still in the 
process of forming their capacitive elements, e.g. by creating pinning centers through crust 



fract ure, while older pulsars have mostly completed the task (lAlpar et al.lll996l : IWong et al. 



200ll ). However, the full picture is more co mplicated. Vela's quasiperiodic avalanches point 
to a richly connected network of reservoirs (lAlpar et al.lll996l ). yet its aggregated value i^gutch 
is relatively low. On the other hand, the other quasiperiodic glitcher, PSR J0537— 6910, 
is relatively young (tc = 4.9 kyr); how did it form a richly connected reservoir network so 
quickly? And, if its network is so richly connected, why is its aggregated i^gutch value so 
low? Likewise, PSR J0358-I-5413 is the oldest object in the sample (tc = 560 kyr), yet its 
e value arguably points to a dearth of capacitive elements, characteristic of a young object. 
There are no obvious grounds (e.g. quasiperiodicity) on which to treat PSR J0358-I-5413 as 
exceptional. 



^^The aggregate value t'gutch (jLvne et al ]l2000l) . binned over semi-decades in z>, effectively averages together 
different pulsars. While this approach reduces the formal error bar on i^gutch, its physical interpretation is 
less straightforward, given the likelihood that e is different in different pulsars. 
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Do all pulsars glitch eventually? It has been speculated in the past that there is some- 
thing special physically about the minority of pulsars that do glitch. While it is impossible 
to reject this hypothesis unequivocally with the data at hand, the results presented here 
suggest that all pulsars are capable of glitching. However, most do so infrequently (low A) 
and hence have not been detected during the last four decades of timing experiments. We 
find that up to ~ 30% of the pulsar population can glitch at rates lower than Amin = 0.25 
yr^^ and still conform with the measured aggregate waiting-time distribution. 

Once verified, the claimed Poissonian nature of the glitch mechanism can be invoked 
to exclude broad classes of glitch theories, e.g. those that rely on 'defects' or 'turbulence' at 
special locations (like the pole), or that involve a pair of dependent events (A. Ma rtin, private 



communication). It is important to interpret aftershocks carefully in this light (IWong et al. 



20011 ). In self-organized critical systems, the excess number of avalanches following a large 
avalanche (over and above the Poissonian baseline following a small avalanche) sc ales in- 
verse ly with the time elapsed, a property known as Omori's law for earthquakes (jJensen 
1998h . 



In this paper, we do not analyze post-glitch relaxation times and glitch-activated changes 
in z> in the context of avalanche proc esses, e.g. the correlation between Az> and the transient 
component of Az/ (jWong et al.ll200ll ). We also assume implicitly that the quantized super- 
fluid vortices in the vortex unpinning model are organized in a rectilinear array, even though 
recent work sugges ts that meridional circula tion destabilizes the array and converts it into 
a turbulent tangle (jPeralta et al.ll2005l . l2006l ). Further study of these matters is deferred to 
future work. 
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Table 1. Parameters of pulsar glitches 
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31.4(2) 


22 








52058(2) 


29(1) 


29 








52802.6(2) 


1.8(7) 


29 


1826- 


-1334 


3 


46507 


2700 


1 








49014 


3060 


1 








53738 


n/a 


29 


1827- 


'0958 


1 


n/a 


n/a 


28 


1833- 


-0559 


1 


52200 


n/a 


29 


1833- 


-0827 


1 


48041 


1864.8 


1 


1835- 


-1106 


1 


52265 


27 


28,29 


1838- 


-0453 


1 


52000 


26 


29 


1841- 


-0425 


1 


53356 


578.5 


29 


1844- 


-0538 


2 


47438 


0.8 


29 
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Table 1 — Continued 



PSR J 






Au/iy 


Ref. 






(MJD) 










17955 


0.5 


29 


1845-0316 


1 


52212.9 


30 


29 


1856+0113 


1 


n/a 


n/a 


28 


1901+0716 


1 


46859 


30 


1 


1902+0615 


4 


48645.11 


0.45 


29 






49441 


0.23 


29 






50311 


0.31 


29 






51165.9 


0.47 


29 


1903+0135 


1 


48634 


n/a 


28,29 


1905-0056 


2 


49385 


n/a 


28,29 






49695.3 


0.21 


28,29 


1908+0909 


2 


52240 


11.8 


29 






53340 


1.7 


29 


1909+0007 


1 


49491.9 


0.72 


29 


1910+0358 


1 


52331 


1.4 


29 


1910-0309 


3 


48241 


0.6 


2 






49219.85 


1.84 


2 






53232.75 


2.66 


29 


1913+0446 


1 


53500 


n/a 


29 


1918+1444 


1 


52285 


2.2 


28,29 


1919+0021 


1 


50174 


1.29 


2 


1922+2018 


1 


n/a 


n/a 


1 


1932+2220 


4 


46900 


4450 


29 






50264 


4457 


2 






52210 


12 


29 






52394 


12 


29 


1946+2611 


1 


53326 


70 


29 


1952+3252 


5 


n/a 


n/a 


28 






51967 


2.25 


31 






52385 


0.72 


31 
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Table 1 — Continued 



T~) O T~) T 


A T 

Ng 


bpocn 


A / 

Au/u 


Kei. 






(MJD) 


(10 








52912 


1.29 


31 






53305 


0.51 


31 


2021+3651 


1 


52630.07 


2587 


18 


2040+1657 


1 


53142 


n/a 


29 


2116+1414 


3 


47972 


0.2 


28,29 






49950 


0.07 


29 






OiOO ( 


n 1 1 


9Q 


2225+6535 


4 


43072 


1707 


1 






51900 


0.14 


31 






52950 


0.08 


31 






53434 


0.19 


31 


2229+6114 


1 


53064 


1133 


29 


2257+5909 


1 


49463.2 


0.92 


2 


2301+5852 


1 


52443.9 


4100 


19 


2330-2005 


1 


n/a 


n/a 


1 


2337+6151 


1 


53639 


20000 


29 
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Table 2. Observing time intervals for pulsars with Ng> 5. 



PSR J 


^min 




^max 






(MJD) 




(MJD) 




0358+5413 


41807 




53546 




0534+2200 


40466 


* 


53476 


* 


0537-6910 


51197 




53952 




0631+1036 


50186 


* 


53621 


* 


0835-4510 


40140 


* 


53960 


* 


1341-6220 


47915 




51022 




1740-3015 


46770 


* 


53190 




1801-2304 


46697 




53356 


* 


1825-0935 


48300 


* 


52803 


* 



Note. — *: Segmented data spans 
are not published, tmm and t^ax are 
estimated by eye from graphed spin- 
down histories, where available, or 
else from the first and last glitches by 
default. 



Table 3: Power-law size distribution parameters for pulsars with A^g > 5. 



PSR J 


a_ 


a 


a+ 


-Pks 


0358+5413 


1.5 


2.4 


5.2 


0.9913 


0534+2200 


1.1 


1.2 


1.4 


0.86 


0537-6910 


0.39 


0.42 


0.43 


0.36 


0631+1036 


1.2 


1.8 


2.7 


0.99896 


0835-4510 


-0.20 


-0.13 


0.18 


0.908 


1341-6220 


1.2 


1.4 


2.1 


0.77 


1740-3015 


0.98 


1.1 


1.3 


0.9920 


1801-2304 


0.092 


0.57 


1.1 


0.99968 


1825-0935 


-0.30 


0.36 


1.0 


0.99904 
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Table 4: Two- component size distribution parameters for quasiperiodic glitcliers. 



PSR J 


a_ 


a 


«+ 






-Pks 


0537-6910 


0.22 


0.44 


0.68 


0.25 ±0.05 


(3.0 ±0.5) X 10-^ 


0.81 


0835-4510 


-0.49 


0.11 


0.44 


0.15 ±0.05 


(2.5 ±0.5) X 10-6 


0.970 



Table 5. Poissonian waiting-time distribution parameters for pulsars with Ng > 5. 



PSR J 


•^min 


•^111111 


Atjiiax 


A_ 


A 


A+ 


-Pks 




(d) 


(d) 


(d) 


(yr-^) 


(yr-^) 


(yr-^) 




0358+5413 


4 


30 


11739 


0.21 


0.57 


1.3 


0.999960 


0534+2200 


2 


18 


13010 


0.57 


0.91 


1.3 


0.982 


0537-6910 


3 


19 


2755 


n/a 


2.6 


n/a 


0.31 


0631+1036 


2 


16 


3435 


0.55 


0.95 


1.9 


0.9970 


0835-4510 


2 


24 


13820 


0.33 


0.35 


0.42 


0.45 


1341-6220 


4 


260 


3107 


1.2 


1.8 


5.6 


0.980 


1740-3015 


4 


100 


6330 


1.2 


1.5 


2.5 


0.928 


1801-2304 


4 


200 


6659 


0.35 


0.55 


0.88 


0.962 


1825-0935 


4 


16 


4503 


0.48 


0.91 


1.8 


0.9989 



Note. — n/a: not applicable, as the best fit yields Pks < 0.32. 



Table 6. Two-component waiting-time distribution parameters for quasiperiodic glitchers. 



PSR J 


A_ 


A 


A+ 




Ate 


Pks 








(yr-^) 




(yr) 




0537-6910 


2.0 


2.6 


3.3 


0.25 ±0.05 


0.3 ±0.1 


0.80 


0835-4510 


0.27 


0.43 


0.62 


0.25 ±0.05 


2.8 ±0.1 


0.968 
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Fig. 1. — Cumulative distribution of fractional glitch sizes Az//z/ in the nine pulsars that 
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Fig. 2. — Cumulative distribution of fractional glitch sizes Au/u for the two pulsars which 
have a quasiperiodic component. The observational data [histogram) are plotted together 
with the best two-component fit {solid curve) given by (jlj), with (Az//z/)inin and (Az//z/)inax 
taken from Table [H and a, Cp and {Au/h')^ taken from Table HI 
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Fig. 3. — Aggregate cumulative distribution of fractional glitch sizes Avjv for all glitch- 
ing pulsars (/e/t panel) and for all glitching pulsars except PSR J0537— 6910 and PSR 
J0835— 4510, which have a quasiperiodic component [right panel). The observational data 
(histogram) are plotted together with the best power-law fit [solid curve) given by ([2]), for 
the best fit parameters in §4.41 
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Cumulative distribution of glitch waiting times At (measured in yr) for the nine 



PSR J0537-6910 



PSR J0835-4510 




time (yr) time (yr) 



Fig. 5. — Cumulative distribution of glitch waiting times At (measured in yr) for the 
two pulsars which have a quasiperiodic component. The observational data {histogram) are 
plotted together with the best two-component fit {solid curve), given by (ITOi) . with Atmin 
taken from Table [T] (twice the epochal uncertainty), Atmax from Table and A, C'^, and At^ 
taken from Table [TJ 
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Fig. 6. — Cumulative distribution of mean glitching rate A (measured in yr~^) for the nine 
pulsars that have glitched more than five times {left panel) and excluding the two quasiperi- 
odic glitchers {right panel), showing the observational data {histogram) and the best fit 
obtained from (fT2|) {solid curve), with (A) = 1.3toljr'^ {left panel) and (A) 
{right panel). 
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Fig. 7. — Aggregate cumulative waiting-time distribution Pagg(At) for all pulsars with Ng > 
2, including {left panel) and excluding {right panel) the two quasiperiodic glitchers. The 
histograms represent the observational data. The dotted curves are obtained by evaluating 
([T^ with ([HD, using (A) = 1.3l|]:^yr-i {left panel) and (A) = 1.21°;^ yr'^ {right panel) 
extracted from Figure O The dashed curves are obtained by evaluating f[T5]) with f[T^ and 
adjusting (A) to maximize Pks when fitting Pagg(At). The solid curves are obtained from 
([ISD, with (A) = 0.54 yr-^ Amm = 0.25 yr^^ {left panel) and (A) = 0.43 yr^^ A^in = 0.25 
{right panel). 
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Fig. 8. — Cumulative distribution of mean glitching rate A (measured in yr"^) for the nine 
pulsars that have glitched more than five times {left panel) and excluding the two quasiperi- 
odic glitchers {right panel), showing the observational data {histogram) and the theoretical 
rate distribution (fill) corresponding to the dashed and solid curves in Figure [71 
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Fig. 9. — Spin-down age (in kyr) versus power-law exponent a for the glitch size distri- 
bution. The error bars indicate the 1-a range of allowable fits according to the K-S test. 
Systematic differences between Tc and true age are not quantified here. Solid (open) tri- 
angles symbolize quasiperiodic glitchers with A^^g > 5, to which we fit a two-component 
(one-component) P(Ai//i/), as in Table H] ([3]). Squares symbolize aperiodic glitchers with 
A^g > 5, to which we fit a power-law P(Az//z/), as in Table [31 
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Fig. 10. — Mean glitching rate A (in yr~^) versus spin-down age (in kyr). The error bars 
indicate the 1-a range of allowable fits according to the K-S test. Systematic differences be- 
tween Tc and true age are not quantified here. Solid (open) triangles symbolize quasiperiodic 
glitchers with Ng > 5, to which we fit a two-component (one-component) P(A, At), as in 
Table [6] ([5]) . Squares symbolize aperiodic glitchers with A'^g > 5, to which we fit a Poissonian 
P(A, At) as in Table E 



